Processing Pipeline for Hematopoesis Data from Popescu et al (2019)

Bacher Group

Author
Affiliation

Jack R. Leary

University of Florida

Published

October 15, 2024

1 Libraries

Before we do anything else we load in the necessary software packages in both R & Python.

1.1 R

Code
library(dplyr)       # data manipulation
library(Seurat)      # scRNAseq tools
library(scLANE)      # trajectory DE
library(ggplot2)     # pretty plots
library(biomaRt)     # gene annotation
library(tradeSeq)    # more trajectory DE
library(patchwork)   # plot alignment
library(slingshot)   # pseudotime estimation
library(reticulate)  # Python interface
select <- dplyr::select
rename <- dplyr::rename

1.2 Python

Code
import scvi                                                     # integration
import warnings                                                 # filter out warnings
import numpy as np                                              # linear algebra tools
import scanpy as sc                                             # scRNA-seq tools
import pandas as pd                                             # DataFrames
import anndata as ad                                            # scRNA-seq data structures
import scvelo as scv                                            # RNA velocity
import cellrank as cr                                           # cell fate estimation
import matplotlib.pyplot as plt                                 # pretty plots
from scipy.io import mmread, mmwrite                            # sparse matrix IO
from cellrank.estimators import GPCCA                           # CellRank estimator
from scipy.sparse import coo_matrix, csr_matrix                 # sparse matrices
from cellrank.kernels import PseudotimeKernel, CytoTRACEKernel  # CellRank kernels

Here we set a few global variables to reduce the amount of printed output caused by our Python code.

Code
warnings.simplefilter('ignore', category=UserWarning)
sc.settings.verbosity = 0
scv.settings.verbosity = 0
cr.settings.verbosity = 0

2 Visualization tools

To make our visualizations prettier we’ll define some consistent color palettes, and set a theme for matplotlib that matches ggplot2::theme_classic().

2.1 Color palettes

Code
palette_heatmap <- paletteer::paletteer_d("MetBrewer::Cassatt1", direction = -1)
palette_celltype <- as.character(paletteer::paletteer_d("ggsci::category10_d3"))[1:6]
names(palette_celltype) <- c("HSC", "Kupffer cell", "Monocyte-macrophage", "Monocyte", "Monocyte precursor", "Neutrophil-myeloid progenitor")
palette_cluster <- paletteer::paletteer_d("ggsci::default_igv")

2.2 Theme for matplotlib

Code
base_size = 12
plt.rcParams.update({
    # font
    'font.size': base_size, 
    'font.weight': 'normal',
    # figure
    'figure.dpi': 320, 
    'figure.edgecolor': 'white', 
    'figure.facecolor': 'white', 
    'figure.figsize': (6, 4), 
    'figure.constrained_layout.use': True,
    # axes
    'axes.edgecolor': 'black',
    'axes.grid': False,
    'axes.labelpad': 2.75,
    'axes.labelsize': base_size * 0.8,
    'axes.linewidth': 1.5,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.titlelocation': 'left',
    'axes.titlepad': 11,
    'axes.titlesize': base_size,
    'axes.titleweight': 'normal',
    'axes.xmargin': 0.1, 
    'axes.ymargin': 0.1, 
    # legend
    'legend.borderaxespad': 1,
    'legend.borderpad': 0.5,
    'legend.columnspacing': 2,
    'legend.fontsize': base_size * 0.8,
    'legend.frameon': False,
    'legend.handleheight': 1,
    'legend.handlelength': 1.2,
    'legend.labelspacing': 1,
    'legend.title_fontsize': base_size, 
    'legend.markerscale': 1.5
})

3 Helper functions

We first define a utility function to make our plot legends cleaner to read & easier to make.

Code
guide_umap <- function(key.size = 4) {
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = key.size,
                                                                    alpha = 1, 
                                                                    stroke = 0.25)))
}

4 Data

4.1 Import hematopoesis data

We start by importing the AnnData object we downloaded from the Human Developmental Cell Atlas portal; it contains a variety of celltypes beyond just the myeloid lineage which we’re interested.

Code
ad_liver = ad.read_h5ad('../../data/hematopoesis/fetal_liver_alladata.h5ad')

We subset to just the myeloid lineage, and remove cells with a percentage mitochondrial reads greater than 15%.

Code
ad_myeloid = ad_liver[ad_liver.obs['cell.labels'].isin(['HSC_MPP', 'Neutrophil-myeloid progenitor', 'Monocyte precursor', 'Monocyte', 'Mono-Mac', 'Kupffer Cell'])]
ad_myeloid.obs['celltype'] = (
    ad_myeloid.obs['cell.labels']
              .map(lambda x: {'HSC_MPP': 'HSC', 'Mono-Mac': 'Monocyte-macrophage', 'Kupffer Cell': 'Kupffer cell'}.get(x, x))
              .astype('category')
)
ad_myeloid = ad_myeloid[ad_myeloid.obs['percent.mito'] < 0.15]
ad_myeloid.layers['counts'] = ad_myeloid.X.copy() 
ad_myeloid.layers['raw_counts'] = ad_myeloid.X.copy() 
ad_myeloid.layers['spliced'] = ad_myeloid.X.copy() 
ad_myeloid.layers['unspliced'] = ad_myeloid.X.copy()
ad_myeloid.raw = ad_myeloid

4.2 Preprocessing

We’ll begin our analysis by running our cells through a standard preprocessing pipeline composed of QC, integration, normalization, HVG selection, dimension reduction, & clustering.

4.2.1 QC

We filter out cells with a sequencing depth of less than 1,000, then remove out genes expressed in less than 5 cells.

Code
sc.pp.filter_cells(ad_myeloid, min_counts=1000)
sc.pp.filter_genes(ad_myeloid, min_cells=5)

4.2.2 HVG selection

Here we select the top 4000 most highly-variable genes (HVGs), which we’ll use as input to PCA, integration, etc.

Code
sc.pp.highly_variable_genes(
    ad_myeloid, 
    n_top_genes=4000, 
    flavor='seurat_v3', 
    layer='counts',
    subset=True
)

4.2.3 Integration with scVI

Using the scVI library, we perform celltype-aware integration across batches by training a variational autoencoder (VAE) for 150 epochs. This provides us with a 20-dimensional integrated latent space embedding that is (hopefully) free of batch effects.

Code
scvi.settings.seed = 312
Code
scvi.settings.num_threads = 16
scvi.model.SCVI.setup_anndata(
    ad_myeloid, 
    layer='counts', 
    batch_key='fetal.ids', 
    labels_key='celltype'
)
int_model = scvi.model.SCVI(
    ad_myeloid, 
    n_layers=2, 
    n_hidden=72, 
    n_latent=20, 
    gene_likelihood='nb', 
    dispersion='gene-label'
)
int_model.train(
    early_stopping=True,
    accelerator='cpu', 
    max_epochs=150
)
Code
ad_myeloid.obsm['X_scVI'] = int_model.get_latent_representation()
Code
sc.pl.embedding(
    ad_myeloid, 
    basis='scVI', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    show=False
)
plt.gca().set_xlabel('scVI 1')
plt.gca().set_ylabel('scVI 2')
plt.show()
Figure 1: scVI embedding colored by celltype.

4.2.4 Normalization

We next depth- and log1p-normalize the mRNA counts matrix.

Code
ad_myeloid.X = sc.pp.normalize_total(ad_myeloid, target_sum=1e4, inplace=False)['X']
sc.pp.log1p(ad_myeloid)

4.2.5 PCA embedding

Using PCA we generate a 50-dimensional linear reduction of the normalized counts.

Code
sc.pp.scale(ad_myeloid)
sc.tl.pca(
    ad_myeloid, 
    n_comps=50, 
    random_state=312, 
    use_highly_variable=True
)

Plotting the first two PCs shows clear separation by celltype, though variation in the immature celltypes seems less clear.

Code
sc.pl.embedding(
    ad_myeloid, 
    basis='pca', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 2: PCA embedding colored by celltype.

4.2.6 SNN graph estimation

Next, we identify the 50 nearest-neighbors (NNs) for each cell using the cosine distance in the scVI latent space.

Code
sc.pp.neighbors(
    ad_myeloid, 
    n_neighbors=50,
    n_pcs=None,  
    metric='cosine', 
    random_state=312, 
    use_rep='X_scVI'
)

4.2.7 Clustering

Using the Leiden algorithm we partition the SNN graph into clusters.

Code
sc.tl.leiden(
    ad_myeloid, 
    resolution=0.3, 
    random_state=312
)

The identified clusters roughly match our celltypes.

Code
sc.pl.embedding(
    ad_myeloid, 
    basis='pca', 
    color='leiden',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 3: PCA embedding colored by Leiden cluster.

4.2.8 Moment-based imputation

Moving on, we create smoothed versions of the counts across each cell’s 30 NNs.

Code
scv.pp.moments(
    ad_myeloid, 
    n_pcs=None,
    use_rep='X_scVI', 
    n_neighbors=30
)

4.2.9 Undirected PAGA embedding

We use the PAGA algorithm to generate a graph abstraction of the relationships between celltypes.

Code
sc.tl.paga(ad_myeloid, groups='celltype')

Visualizing the results shows a pretty linear structure that matches what we expect biologically - stem cells at the root, monocytes & macrophages in the middle, and Kupffer cells at the end.

Code
sc.pl.paga(
    ad_myeloid, 
    fontoutline=True, 
    fontsize=12, 
    frameon=True, 
    random_state=312, 
    show=False
)
plt.gca().set_xlabel('PAGA 1')
plt.gca().set_ylabel('PAGA 2')
plt.show()
Figure 4: PAGA embedding colored by celltype. Lines are weighted by (undirected) connectivity strength.

4.2.10 UMAP embedding

We then create a nonlinear dimension reduction of the data using the UMAP algorithm.

Code
sc.tl.umap(ad_myeloid, random_state=312)

In UMAP space the transitions between celltypes are clear and well-represented.

Code
sc.pl.embedding(
    ad_myeloid, 
    basis='umap', 
    color='celltype',
    title='',
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 5: UMAP embedding colored by celltype.

4.2.11 Force-directed graph embedding

Using the ForceAtlas2 algorithm we generate a 2-dimensional graph layout of our cells.

Code
sc.tl.draw_graph(
    ad_myeloid, 
    layout='fa', 
    random_state=312
)

The force-directed graph embedding is more visually appealing and cleaner than the UMAP embedding, so we’ll most likely use it instead going forwards.

Code
sc.pl.draw_graph(
    ad_myeloid, 
    color='celltype', 
    title='', 
    size=30, 
    alpha=0.5, 
    show=False, 
    frameon=True
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 6: Force-directed graph embedding colored by celltype.

4.2.12 Diffusion map embedding

Our last embedding will be generated using diffusion maps, which explicitly attempts to preserve transitional structures in the data.

Code
sc.tl.diffmap(
    ad_myeloid, 
    random_state=312, 
    n_comps=16
)
ad_myeloid.obsm['X_diffmap_old'] = ad_myeloid.obsm['X_diffmap']
ad_myeloid.obsm['X_diffmap'] = ad_myeloid.obsm['X_diffmap'][:, 1:] 

While the diffusion map embedding captures some characteristics of the data, it masks variation in the mature macrophage celltypes.

Code
sc.pl.embedding(
    ad_myeloid, 
    basis='diffmap', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 7: Diffusion map embedding colored by celltype.

4.2.13 Diffusion pseudotime estimation

We next compute a diffusion pseudotime estimate for each cell using the first diffusion component.

Code
ad_myeloid.uns['iroot'] = np.argmin(ad_myeloid.obsm['X_diffmap'][:, 0])
sc.tl.dpt(ad_myeloid, n_dcs=1)

Overall the diffusion pseudotime captures myeloid differentiation progression well, but it assigns very similar values to all the mature Kupffer cells which isn’t ideal.

Code
sc.pl.embedding(
    ad_myeloid, 
    basis='diffmap', 
    color='dpt_pseudotime',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 8: Diffusion map embedding colored by diffusion pseudotime.

4.3 Save preprocessed data

We save the processed AnnData object to disk.

Code
ad_myeloid.write_h5ad('../../data/hematopoesis/ad_myeloid.h5ad')

We also save the embeddings, counts matrix, metadata, etc. so that we can create a Seurat object.

Code
ad_myeloid.obs.reset_index(inplace=False).rename(columns={'index': 'cell'}).to_csv('../../data/hematopoesis/conversion_files/cell_metadata.csv')
ad_myeloid.var.reset_index(inplace=False).rename(columns={'index': 'gene'}).to_csv('../../data/hematopoesis/conversion_files/gene_metadata.csv')
pd.DataFrame(ad_myeloid.obsm['X_draw_graph_fa']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/graph_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_umap']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/umap_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_pca']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/pca_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_scVI']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/scvi_embed.csv')
pd.DataFrame(ad_myeloid.obsm['X_diffmap']).reset_index(inplace=False).to_csv('../../data/hematopoesis/conversion_files/diffmap_embed.csv')
mmwrite('../../data/hematopoesis/conversion_files/raw_counts.mtx', a=ad_myeloid.layers['raw_counts'])

5 Analysis

5.1 Read data into R

We start our analysis by importing the preprocessed data from Python into R. Expand the code block for details.

Code
cell_metadata <- readr::read_csv("../../data/hematopoesis/conversion_files/cell_metadata.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1) %>% 
                 as.data.frame() %>% 
                 magrittr::set_rownames(.$cell)
gene_metadata <- readr::read_csv("../../data/hematopoesis/conversion_files/gene_metadata.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1)
embed_PCA <- readr::read_csv("../../data/hematopoesis/conversion_files/pca_embed.csv",
                             show_col_types = FALSE, 
                             col_select = -c(1:2)) %>% 
             as.data.frame() %>% 
             magrittr::set_colnames(paste0("PC", 1:ncol(.))) %>% 
             magrittr::set_rownames(cell_metadata$cell)
embed_PCA <- CreateDimReducObject(embeddings = as.matrix(embed_PCA),
                                  key = "PC_", 
                                  global = TRUE, 
                                  assay = "RNA")
embed_FA <- readr::read_csv("../../data/hematopoesis/conversion_files/graph_embed.csv", 
                            show_col_types = FALSE, 
                            col_select = -c(1:2)) %>% 
            as.data.frame() %>% 
            magrittr::set_colnames(paste0("FA", 1:ncol(.))) %>% 
            magrittr::set_rownames(cell_metadata$cell)
embed_FA <- CreateDimReducObject(embeddings = as.matrix(embed_FA),
                                 key = "FA_", 
                                 global = TRUE, 
                                 assay = "RNA")
embed_UMAP <- readr::read_csv("../../data/hematopoesis/conversion_files/umap_embed.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("UMAP", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_UMAP <- CreateDimReducObject(embeddings = as.matrix(embed_UMAP),
                                   key = "UMAP_", 
                                   global = TRUE, 
                                   assay = "RNA")
embed_diffmap <- readr::read_csv("../../data/hematopoesis/conversion_files/diffmap_embed.csv", 
                                 show_col_types = FALSE, 
                                  col_select = -c(1:2)) %>% 
                 as.data.frame() %>% 
                 magrittr::set_colnames(paste0("DC", 1:ncol(.))) %>% 
                 magrittr::set_rownames(cell_metadata$cell)
embed_diffmap <- CreateDimReducObject(embeddings = as.matrix(embed_diffmap),
                                      key = "DC_", 
                                      global = TRUE, 
                                      assay = "RNA")
embed_scVI <- readr::read_csv("../../data/hematopoesis/conversion_files/scvi_embed.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("scVI", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_scVI <- CreateDimReducObject(embeddings = as.matrix(embed_scVI),
                                   key = "scVI_", 
                                   global = TRUE, 
                                   assay = "RNA")
raw_counts <- Matrix::t(Matrix::readMM("../../data/hematopoesis/conversion_files/raw_counts.mtx"))
colnames(raw_counts) <- cell_metadata$cell
rownames(raw_counts) <- gene_metadata$gene
raw_counts <- CreateAssayObject(counts = raw_counts,
                                min.cells = 0,
                                min.features = 0, 
                                key = "rna")

We can now create a Seurat object, to which we add our various embeddings and metadata. After adding all the various bits and pieces, we normalize the raw counts and re-identify HVGs.

Code
seu_myeloid <- CreateSeuratObject(raw_counts, 
                                  assay = "RNA",
                                  meta.data = cell_metadata, 
                                  project = "hematopoesis", 
                                  min.cells = 0, 
                                  min.features = 0)
seu_myeloid@meta.data <- mutate(seu_myeloid@meta.data, 
                                celltype = factor(celltype, levels = c("HSC", 
                                                                       "Kupffer cell", 
                                                                       "Monocyte-macrophage", 
                                                                       "Monocyte", 
                                                                       "Monocyte precursor", 
                                                                       "Neutrophil-myeloid progenitor")))
seu_myeloid@reductions$pca <- embed_PCA
seu_myeloid@reductions$diffmap <- embed_diffmap
seu_myeloid@reductions$fa <- embed_FA
seu_myeloid@reductions$umap <- embed_UMAP
seu_myeloid@reductions$scvi <- embed_scVI
seu_myeloid <- NormalizeData(seu_myeloid, verbose = FALSE) %>% 
               FindVariableFeatures(selection.method = "vst", verbose = FALSE)

5.2 Pseudotime estimation

Using Slingshot we identify a single lineage of pseudotime across the force-directed graph embedding.

Code
sling_res <- slingshot(seu_myeloid@reductions$fa@cell.embeddings[, 1:2], 
                       clusterLabels = as.factor(seu_myeloid$celltype), 
                       start.clus = c("HSC"), 
                       approx_points = 1000)
sling_curves <- slingCurves(sling_res, as.df = TRUE)
sling_mst <- slingMST(sling_res, as.df = TRUE)
sling_pt <- slingPseudotime(sling_res) %>%
            as.data.frame() %>%
            mutate(cell = rownames(.), .before = 1) %>%
            rename(PT = Lineage1) %>% 
            mutate(PT = (PT - min(PT)) / (max(PT) - min(PT)))

Visualizing the results from Slingshot shows us that our biological expectations have been met.

Code
p0 <- data.frame(Embeddings(seu_myeloid, "fa")) %>% 
      mutate(celltype = seu_myeloid$celltype) %>% 
      ggplot(aes(x = FA_1, y = FA_2, color = celltype)) + 
      geom_point(size = 1.5, 
                 alpha = 0.5, 
                 stroke = 0) + 
      geom_path(data = sling_mst, mapping = aes(x = FA_1, y = FA_2, group = Lineage), 
                linewidth = 1.25, 
                color = "black") + 
      geom_point(data = sling_mst, mapping = aes(x = FA_1, y = FA_2, fill = Cluster), 
                color = "black", 
                shape = 21, 
                size = 4.5, 
                stroke = 1.25, 
                show.legend = FALSE) + 
      scale_color_manual(values = palette_celltype) + 
      scale_fill_manual(values = palette_celltype) + 
      labs(x = "FA 1", y = "FA 2") +
      theme_scLANE(umap = TRUE) +  
      theme(legend.title = element_blank()) + 
      guide_umap()
p1 <- data.frame(Embeddings(seu_myeloid, "fa")) %>% 
      mutate(celltype = seu_myeloid$celltype) %>% 
      ggplot(aes(x = FA_1, y = FA_2, color = celltype)) + 
      geom_point(size = 1.5, 
                 alpha = 0.5, 
                 stroke = 0) + 
      geom_path(data = sling_curves,
                mapping = aes(x = FA_1, y = FA_2, group = Lineage), 
                color = "black", 
                linewidth = 1.5, 
                alpha = 0.5, 
                lineend = "round") + 
      scale_color_manual(values = palette_celltype) + 
      labs(x = "FA 1", y = "FA 2") + 
      theme_scLANE(umap = TRUE) + 
      theme(legend.title = element_blank()) + 
      guide_umap()
p2 <- (p0 / p1) + plot_layout(guides = "collect", axes = "collect")
p2
Figure 9: Force-directed graph embedding overlaid with the MST and principal curves estimated by Slingshot.

5.3 Trajectory DE testing with scLANE

We use the scLANE package to perform trajectory DE testing across pseudotime time for the top 4000 HVGs.

Code
cell_offset <- createCellOffset(seu_myeloid)
pt_df <- data.frame(PT = sling_pt$PT)
candidate_genes <- chooseCandidateGenes(seu_myeloid, 
                                        group.by.subject = FALSE, 
                                        n.desired.genes = 4000L)
Warning in asMethod(object): sparse->dense coercion: allocating vector of size
1.1 GiB
Code
scLANE_models <- testDynamic(seu_myeloid, 
                             pt = pt_df, 
                             genes = candidate_genes, 
                             size.factor.offset = cell_offset, 
                             n.cores = 16L,
                             verbose = FALSE)
Warning in asMethod(object): sparse->dense coercion: allocating vector of size
1.1 GiB
scLANE testing in GLM mode completed for 4000 genes across 1 lineage in 10.232 hours
Code
scLANE_de_res <- getResultsDE(scLANE_models)

5.4 tradeSeq DE analysis

Code
ts_start <- Sys.time()
bioc_par <- BiocParallel::MulticoreParam(workers = 16L, RNGseed = 312)
RNA_counts <- as.matrix(seu_myeloid@assays$RNA$counts)[candidate_genes, ]
Warning in asMethod(object): sparse->dense coercion: allocating vector of size
1.1 GiB
Code
k_eval <- evaluateK(RNA_counts, 
                    pseudotime = pt_df, 
                    cellWeights = matrix(rep(1, nrow(pt_df)), ncol = 1), 
                    offset = log(1 / cell_offset), 
                    k = 3:10, 
                    plot = FALSE, 
                    nGenes = 500, 
                    verbose = FALSE, 
                    parallel = TRUE, 
                    BPPARAM = bioc_par)
best_k <- c(3:10)[which.min(abs(colMeans(k_eval - rowMeans(k_eval))))]  # choose k w/ lowest MAD from mean AIC 
ts_models <- fitGAM(RNA_counts, 
                    pseudotime = pt_df, 
                    cellWeights = matrix(rep(1, nrow(pt_df)), ncol = 1), 
                    offset = log(1 / cell_offset), 
                    nknots = best_k, 
                    sce = FALSE, 
                    parallel = TRUE, 
                    verbose = FALSE, 
                    BPPARAM = bioc_par)
BiocParallel::bpstop(bioc_par)
ts_de_res <- associationTest(ts_models, global = TRUE) %>% 
                             arrange(desc(waldStat)) %>% 
                             mutate(gene = rownames(.), 
                                    pvalue_adj = p.adjust(pvalue, method = "fdr"), 
                                    gene_dynamic_overall = if_else(pvalue_adj < 0.01, 1, 0)) %>% 
                             relocate(gene)
ts_end <- Sys.time()
ts_diff <- ts_end - ts_start

The runtime for tradeSeq with 5 knots per gene is:

Code
ts_diff
Time difference of 2.700508 hours

5.5 Cell fate analysis

Lastly, we’ll perform an analysis of cell fate specification using the CellRank Python package.

5.5.1 CytoTRACE kernel

We begin by computing a CytoTRACE kernel; CytoTRACE uses the number of genes expressed in a cell as a proxy measurement for differentiation potential, with the assumption being that cells expressing more genes are more likely to be progenitor celltypes & vice versa.

Code
ctk = CytoTRACEKernel(ad_myeloid).compute_cytotrace()

We plot the CytoTRACE scores below. Interestingly, in addition to the expected result of HSCs having high CytoTRACE scores we see a subcluster of Kupffer cells with high scores as well.

Code
sc.pl.embedding(
    ad_myeloid,
    basis='draw_graph_fa',
    color=['ct_score', 'ct_pseudotime'],
    color_map='magma', 
    show=False, 
    size=30, 
    alpha=0.5, 
    title=['CytoTRACE score', 'CytoTRACE Pseudotime']
)
plt.show()
Figure 10: Force-directed graph embedding colored by CytoTRACE score and pseudotime. Higher scores indicate higher differentiation potential.

Next, we estimate a cell-cell transition probability matrix.

Code
ctk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the projection of that matrix on our force-directed graph embedding shows messy / reversed directionality. This indicates that CytoTRACE might not work very well on this dataset, and at the very least should be augmented with another data source.

Code
ctk.plot_projection(
  basis='draw_graph_fa', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.5, 
  linewidth=2, 
  frameon=True, 
  show=False
)
/blue/rbacher/j.leary/py_envs/scLANE_env2/lib/python3.10/site-packages/scvelo/plotting/utils.py:63: DeprecationWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, pd.CategoricalDtype) instead
  return isinstance(c, str) and c in data.obs.keys() and cat(data.obs[c])
/blue/rbacher/j.leary/py_envs/scLANE_env2/lib/python3.10/site-packages/scvelo/plotting/utils.py:63: DeprecationWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, pd.CategoricalDtype) instead
  return isinstance(c, str) and c in data.obs.keys() and cat(data.obs[c])
Code
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()  
Figure 11: Force-directed graph embedding streamline plot showing differentiation directions as inferred from CytoTRACE scores.

5.5.2 Pseudotime kernel

Next, we use the pseudotime estimates from Slingshot to create a pseudotime-based kernel and estimate another cell-cell transition probability matrix.

Code
ad_myeloid.obs['sling_PT'] = r.sling_pt['PT'].tolist()
pk = PseudotimeKernel(ad_myeloid, time_key='sling_PT')
pk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the projection of the pseudotime kernel’s transition probability matrix shows a much more reasonable portrayal of differentiation directionality.

Code
pk.plot_projection(
  basis='draw_graph_fa', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 12: Force-directed graph embedding overlaid with streamline arrows derived from Slingshot pseudotime.

5.5.3 Combined kernel

We combine the two kernel using unequal weighting, with the pseudotime kernel being given priority due to its better overall capture of differentiation.

Code
ck = 0.2 * ctk + 0.8 * pk

We visualize the combined kernel’s projection below. It seems to reliably portray the direction of differentiation in our dataset.

Code
ck.plot_projection(
  basis='draw_graph_fa', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 13: Force-directed graph embedding overlaid with streamline arrows derived from the combined kernel.

We can also simulate random walks on the high-dimensional cell-cell graph starting from the HSC population. We see that nearly all of the random walks terminate in the Kupffer cell population.

Code
ck.plot_random_walks(
    n_sims=200,
    start_ixs={'celltype': 'HSC'},
    basis='draw_graph_fa',
    color='celltype',
    legend_loc='right margin',
    seed=312, 
    frameon=True, 
    title='', 
    linewidth=0.5, 
    linealpha=0.25, 
    size=30, 
    alpha=0.5, 
    ixs_legend_loc='lower', 
    show_progress_bar=False
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 14: Force-directed graph embedding overlaid with random walks along the cell-cell graph derived from the combined kernel. Starting points are highlighted in black and ending points in yellow.

6 Save data

Finally, we save our Seurat object, the output from Slingshot, and the models from scLANE and tradeSeq.

Code
saveRDS(seu_myeloid, file = "../../data/hematopoesis/seu_myeloid.Rds")
saveRDS(sling_res, file = "../../data/hematopoesis/sling_res.Rds")
saveRDS(scLANE_models, file = "../../data/hematopoesis/scLANE_models.Rds")
saveRDS(ts_models, file = "../../data/hematopoesis/ts_models.Rds")

7 Session info

Code
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.3 (2024-02-29)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/New_York
 date     2024-10-15
 pandoc   2.9.2.1 @ /usr/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date (UTC) lib source
 abind                  1.4-8      2024-09-12 [2] CRAN (R 4.3.3)
 AnnotationDbi          1.62.2     2023-07-02 [2] Bioconductor
 backports              1.5.0      2024-05-23 [1] CRAN (R 4.2.3)
 bigassertr             0.1.6      2023-01-10 [1] CRAN (R 4.3.2)
 bigparallelr           0.3.2      2021-10-02 [1] CRAN (R 4.3.2)
 bigstatsr              1.5.12     2022-10-14 [1] CRAN (R 4.3.2)
 Biobase              * 2.60.0     2023-04-25 [2] Bioconductor
 BiocFileCache          2.13.1     2024-10-05 [1] Github (Bioconductor/BiocFileCache@30abbe4)
 BiocGenerics         * 0.46.0     2023-04-25 [2] Bioconductor
 BiocParallel           1.34.2     2023-05-22 [2] Bioconductor
 biomaRt              * 2.56.1     2023-06-09 [1] Bioconductor
 Biostrings             2.68.1     2023-05-16 [2] Bioconductor
 bit                    4.5.0      2024-09-20 [2] CRAN (R 4.3.3)
 bit64                  4.5.2      2024-09-22 [2] CRAN (R 4.3.3)
 bitops                 1.0-8      2024-07-29 [2] CRAN (R 4.3.3)
 blob                   1.2.4      2023-03-17 [2] CRAN (R 4.3.1)
 boot                   1.3-31     2024-08-28 [4] CRAN (R 4.3.3)
 broom                  1.0.6      2024-05-17 [1] CRAN (R 4.2.3)
 broom.mixed            0.2.9.5    2024-04-01 [1] CRAN (R 4.2.3)
 cachem                 1.1.0      2024-05-16 [2] CRAN (R 4.3.3)
 cli                    3.6.2      2023-12-11 [1] CRAN (R 4.2.3)
 cluster                2.1.6      2023-12-01 [4] CRAN (R 4.3.2)
 coda                   0.19-4.1   2024-01-31 [2] CRAN (R 4.3.2)
 codetools              0.2-20     2024-03-31 [4] CRAN (R 4.3.3)
 colorspace             2.1-1      2024-07-26 [2] CRAN (R 4.3.3)
 cowplot                1.1.3      2024-01-22 [2] CRAN (R 4.3.3)
 crayon                 1.5.3      2024-06-20 [2] CRAN (R 4.3.3)
 curl                   5.2.3      2024-09-20 [2] CRAN (R 4.3.3)
 data.table             1.16.0     2024-08-27 [2] CRAN (R 4.3.3)
 DBI                    1.2.3      2024-06-02 [2] CRAN (R 4.3.3)
 dbplyr                 2.5.0      2024-03-19 [1] CRAN (R 4.3.3)
 DelayedArray           0.26.7     2023-07-28 [2] Bioconductor
 DelayedMatrixStats     1.22.6     2023-08-28 [2] Bioconductor
 deldir                 2.0-4      2024-02-28 [2] CRAN (R 4.3.3)
 digest                 0.6.35     2024-03-11 [1] CRAN (R 4.2.3)
 doParallel             1.0.17     2022-02-07 [1] CRAN (R 4.2.3)
 doSNOW                 1.0.20     2022-02-04 [1] CRAN (R 4.3.2)
 dotCall64              1.1-1      2023-11-28 [2] CRAN (R 4.3.1)
 dplyr                * 1.1.4      2023-11-17 [1] CRAN (R 4.2.3)
 edgeR                  3.42.4     2023-05-31 [2] Bioconductor
 emmeans                1.10.2     2024-05-20 [1] CRAN (R 4.3.3)
 estimability           1.5.1      2024-05-12 [1] CRAN (R 4.3.3)
 evaluate               0.24.0     2024-06-10 [1] CRAN (R 4.3.3)
 fansi                  1.0.6      2023-12-08 [2] CRAN (R 4.3.1)
 farver                 2.1.2      2024-05-13 [1] CRAN (R 4.2.3)
 fastDummies            1.7.4      2024-08-16 [2] CRAN (R 4.3.3)
 fastmap                1.2.0      2024-05-15 [2] CRAN (R 4.3.3)
 ff                     4.0.12     2024-01-12 [1] CRAN (R 4.3.2)
 filelock               1.0.3      2023-12-11 [2] CRAN (R 4.3.3)
 fitdistrplus           1.2-1      2024-07-12 [2] CRAN (R 4.3.3)
 flock                  0.7        2016-11-12 [1] CRAN (R 4.3.2)
 forcats                1.0.0      2023-01-29 [1] CRAN (R 4.3.2)
 foreach                1.5.2      2022-02-02 [1] CRAN (R 4.2.3)
 furrr                  0.3.1      2022-08-15 [1] CRAN (R 4.3.2)
 future                 1.33.2     2024-03-26 [1] CRAN (R 4.2.3)
 future.apply           1.11.2     2024-03-28 [2] CRAN (R 4.3.3)
 gamlss                 5.4-22     2024-03-20 [1] CRAN (R 4.3.2)
 gamlss.data            6.0-6      2024-03-14 [1] CRAN (R 4.3.2)
 gamlss.dist            6.1-1      2023-08-23 [1] CRAN (R 4.3.2)
 geeM                   0.10.1     2018-06-18 [1] CRAN (R 4.3.2)
 generics               0.1.3      2022-07-05 [1] CRAN (R 4.2.0)
 GenomeInfoDb         * 1.36.4     2023-10-02 [2] Bioconductor
 GenomeInfoDbData       1.2.10     2023-07-26 [2] Bioconductor
 GenomicRanges        * 1.52.1     2023-10-08 [2] Bioconductor
 ggplot2              * 3.5.1      2024-04-23 [1] CRAN (R 4.2.3)
 ggrepel                0.9.2      2022-11-06 [1] CRAN (R 4.2.2)
 ggridges               0.5.6      2024-01-23 [2] CRAN (R 4.3.3)
 glm2                 * 1.2.1      2018-08-11 [1] CRAN (R 4.3.2)
 glmmTMB                1.1.9      2024-03-20 [1] CRAN (R 4.3.3)
 globals                0.16.3     2024-03-08 [2] CRAN (R 4.3.3)
 glue                   1.8.0      2024-09-30 [2] CRAN (R 4.3.3)
 goftest                1.2-3      2021-10-07 [2] CRAN (R 4.3.1)
 gridExtra              2.3        2017-09-09 [2] CRAN (R 4.3.1)
 gtable                 0.3.5      2024-04-22 [2] CRAN (R 4.3.3)
 hms                    1.1.3      2023-03-21 [2] CRAN (R 4.3.1)
 htmltools              0.5.8.1    2024-04-04 [1] CRAN (R 4.2.3)
 htmlwidgets            1.6.4      2023-12-06 [2] CRAN (R 4.3.3)
 httpuv                 1.6.15     2024-03-26 [2] CRAN (R 4.3.3)
 httr                   1.4.7      2023-08-15 [2] CRAN (R 4.3.3)
 ica                    1.0-3      2022-07-08 [2] CRAN (R 4.3.1)
 igraph                 2.0.3      2024-03-13 [1] CRAN (R 4.2.3)
 IRanges              * 2.34.1     2023-06-22 [2] Bioconductor
 irlba                  2.3.5.1    2022-10-03 [2] CRAN (R 4.3.1)
 iterators              1.0.14     2022-02-05 [2] CRAN (R 4.3.1)
 jsonlite               1.8.8      2023-12-04 [1] CRAN (R 4.2.3)
 KEGGREST               1.40.1     2023-09-29 [2] Bioconductor
 KernSmooth             2.23-24    2024-05-17 [4] CRAN (R 4.3.3)
 knitr                  1.48       2024-07-07 [2] CRAN (R 4.3.3)
 labeling               0.4.3      2023-08-29 [2] CRAN (R 4.3.1)
 later                  1.3.2      2023-12-06 [2] CRAN (R 4.3.3)
 lattice                0.22-6     2024-03-20 [4] CRAN (R 4.3.3)
 lazyeval               0.2.2      2019-03-15 [2] CRAN (R 4.3.1)
 leiden                 0.4.3.1    2023-11-17 [2] CRAN (R 4.3.1)
 lifecycle              1.0.4      2023-11-07 [2] CRAN (R 4.3.1)
 limma                  3.56.2     2023-06-04 [2] Bioconductor
 listenv                0.9.1      2024-01-29 [2] CRAN (R 4.3.3)
 lme4                   1.1-35.5   2024-07-03 [2] CRAN (R 4.3.3)
 lmtest                 0.9-40     2022-03-21 [2] CRAN (R 4.3.1)
 locfit                 1.5-9.10   2024-06-24 [2] CRAN (R 4.3.3)
 magrittr             * 2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
 MASS                   7.3-60     2023-05-04 [4] CRAN (R 4.3.1)
 Matrix                 1.6-5      2024-01-11 [1] CRAN (R 4.2.3)
 MatrixGenerics       * 1.12.3     2023-07-30 [2] Bioconductor
 matrixStats          * 1.4.1      2024-09-08 [2] CRAN (R 4.3.3)
 memoise                2.0.1      2021-11-26 [2] CRAN (R 4.3.1)
 mgcv                   1.9-1      2023-12-21 [4] CRAN (R 4.3.2)
 mime                   0.12       2021-09-28 [2] CRAN (R 4.3.1)
 miniUI                 0.1.1.1    2018-05-18 [2] CRAN (R 4.3.1)
 minqa                  1.2.7      2024-05-20 [1] CRAN (R 4.2.3)
 munsell                0.5.1      2024-04-01 [2] CRAN (R 4.3.3)
 mvtnorm                1.3-1      2024-09-03 [2] CRAN (R 4.3.3)
 nlme                   3.1-166    2024-08-14 [4] CRAN (R 4.3.3)
 nloptr                 2.1.1      2024-06-25 [2] CRAN (R 4.3.3)
 numDeriv               2016.8-1.1 2019-06-06 [2] CRAN (R 4.3.1)
 paletteer              1.6.0      2024-01-21 [1] CRAN (R 4.2.3)
 parallelly             1.37.1     2024-02-29 [1] CRAN (R 4.2.3)
 patchwork            * 1.2.0      2024-01-08 [1] CRAN (R 4.2.3)
 pbapply                1.7-2      2023-06-27 [2] CRAN (R 4.3.1)
 pillar                 1.9.0      2023-03-22 [2] CRAN (R 4.3.1)
 pkgconfig              2.0.3      2019-09-22 [2] CRAN (R 4.3.1)
 plotly                 4.10.4     2024-01-13 [2] CRAN (R 4.3.3)
 plyr                   1.8.9      2023-10-02 [2] CRAN (R 4.3.3)
 png                    0.1-8      2022-11-29 [2] CRAN (R 4.3.1)
 polyclip               1.10-7     2024-07-23 [2] CRAN (R 4.3.3)
 prettyunits            1.2.0      2023-09-24 [2] CRAN (R 4.3.3)
 princurve            * 2.1.6      2021-01-18 [1] CRAN (R 4.3.2)
 prismatic              1.1.2      2024-04-10 [1] CRAN (R 4.2.3)
 progress               1.2.3      2023-12-06 [2] CRAN (R 4.3.3)
 progressr              0.14.0     2023-08-10 [2] CRAN (R 4.3.2)
 promises               1.3.0      2024-04-05 [2] CRAN (R 4.3.3)
 ps                     1.8.0      2024-09-12 [2] CRAN (R 4.3.3)
 purrr                  1.0.2      2023-08-10 [1] CRAN (R 4.2.3)
 R6                     2.5.1      2021-08-19 [2] CRAN (R 4.3.1)
 RANN                   2.6.2      2024-08-25 [2] CRAN (R 4.3.3)
 rappdirs               0.3.3      2021-01-31 [2] CRAN (R 4.3.1)
 RColorBrewer           1.1-3      2022-04-03 [2] CRAN (R 4.3.1)
 Rcpp                   1.0.13     2024-07-17 [2] CRAN (R 4.3.3)
 RcppAnnoy              0.0.22     2024-01-23 [2] CRAN (R 4.3.3)
 RcppEigen              0.3.4.0.0  2024-02-28 [1] CRAN (R 4.2.3)
 RcppHNSW               0.6.0      2024-02-04 [2] CRAN (R 4.3.3)
 RCurl                  1.98-1.16  2024-07-11 [2] CRAN (R 4.3.3)
 readr                  2.1.5      2024-01-10 [2] CRAN (R 4.3.3)
 rematch2               2.1.2      2020-05-01 [2] CRAN (R 4.3.1)
 reshape2               1.4.4      2020-04-09 [2] CRAN (R 4.3.1)
 reticulate           * 1.35.0     2024-01-31 [1] CRAN (R 4.2.3)
 rlang                  1.1.4      2024-06-04 [1] CRAN (R 4.3.3)
 rmarkdown              2.27       2024-05-17 [1] CRAN (R 4.3.3)
 rmio                   0.4.0      2022-02-17 [1] CRAN (R 4.3.2)
 ROCR                   1.0-11     2020-05-02 [2] CRAN (R 4.3.1)
 RSpectra               0.16-2     2024-07-18 [2] CRAN (R 4.3.3)
 RSQLite                2.3.7      2024-05-27 [2] CRAN (R 4.3.3)
 Rtsne                  0.17       2023-12-07 [2] CRAN (R 4.3.3)
 S4Arrays               1.2.1      2024-03-04 [1] Bioconductor 3.18 (R 4.3.2)
 S4Vectors            * 0.38.2     2023-09-22 [2] Bioconductor
 scales                 1.3.0      2023-11-28 [1] CRAN (R 4.2.3)
 scattermore            1.2        2023-06-12 [2] CRAN (R 4.3.1)
 scLANE               * 0.8.4      2024-10-13 [1] Github (jr-leary7/scLANE@5efbe2c)
 sctransform            0.4.1      2023-10-19 [1] CRAN (R 4.2.3)
 sessioninfo            1.2.2      2021-12-06 [2] CRAN (R 4.3.1)
 Seurat               * 5.1.0      2024-05-10 [1] CRAN (R 4.3.3)
 SeuratObject         * 5.0.2      2024-05-08 [1] CRAN (R 4.3.3)
 shiny                  1.9.1      2024-08-01 [2] CRAN (R 4.3.3)
 SingleCellExperiment * 1.22.0     2023-04-25 [2] Bioconductor
 slingshot            * 2.6.0      2022-11-01 [1] Bioconductor
 snow                   0.4-4      2021-10-27 [2] CRAN (R 4.3.1)
 sp                   * 2.1-4      2024-04-30 [3] CRAN (R 4.3.3)
 spam                   2.10-0     2023-10-23 [2] CRAN (R 4.3.1)
 sparseMatrixStats      1.12.2     2023-07-02 [2] Bioconductor
 spatstat.data          3.1-2      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.explore       3.3-2      2024-08-21 [2] CRAN (R 4.3.3)
 spatstat.geom          3.3-3      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.random        3.3-2      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.sparse        3.1-0      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.univar        3.0-1      2024-09-05 [2] CRAN (R 4.3.3)
 spatstat.utils         3.1-0      2024-08-17 [2] CRAN (R 4.3.3)
 stringi                1.8.4      2024-05-06 [2] CRAN (R 4.3.3)
 stringr                1.5.1      2023-11-14 [2] CRAN (R 4.3.2)
 SummarizedExperiment * 1.30.2     2023-06-06 [2] Bioconductor
 survival               3.7-0      2024-06-05 [2] CRAN (R 4.3.3)
 tensor                 1.5        2012-05-05 [2] CRAN (R 4.3.1)
 tibble                 3.2.1      2023-03-20 [2] CRAN (R 4.3.1)
 tidyr                  1.3.1      2024-01-24 [1] CRAN (R 4.2.3)
 tidyselect             1.2.1      2024-03-11 [2] CRAN (R 4.3.3)
 TMB                    1.9.11     2024-04-03 [1] CRAN (R 4.3.3)
 tradeSeq             * 1.14.0     2023-04-25 [1] Bioconductor
 TrajectoryUtils      * 1.6.0      2022-11-01 [1] Bioconductor
 tzdb                   0.4.0      2023-05-12 [2] CRAN (R 4.3.1)
 utf8                   1.2.4      2023-10-22 [2] CRAN (R 4.3.1)
 uwot                   0.2.2      2024-04-21 [2] CRAN (R 4.3.3)
 vctrs                  0.6.5      2023-12-01 [2] CRAN (R 4.3.1)
 viridis                0.6.5      2024-01-29 [2] CRAN (R 4.3.3)
 viridisLite            0.4.2      2023-05-02 [2] CRAN (R 4.3.1)
 vroom                  1.6.5      2023-12-05 [2] CRAN (R 4.3.3)
 withr                  3.0.0      2024-01-16 [1] CRAN (R 4.2.3)
 xfun                   0.44       2024-05-15 [1] CRAN (R 4.3.3)
 XML                    3.99-0.17  2024-06-25 [2] CRAN (R 4.3.3)
 xml2                   1.3.6      2023-12-04 [2] CRAN (R 4.3.3)
 xtable                 1.8-4      2019-04-21 [2] CRAN (R 4.3.1)
 XVector                0.40.0     2023-04-25 [2] Bioconductor
 yaml                   2.3.10     2024-07-26 [2] CRAN (R 4.3.3)
 zlibbioc               1.46.0     2023-04-25 [2] Bioconductor
 zoo                    1.8-12     2023-04-13 [2] CRAN (R 4.3.1)

 [1] /home/j.leary/r_packages_default
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /blue/rbacher/j.leary/py_envs/scLANE_env2/bin/python
 libpython:      /blue/rbacher/j.leary/py_envs/scLANE_env2/lib/libpython3.10.so
 pythonhome:     /blue/rbacher/j.leary/py_envs/scLANE_env2:/blue/rbacher/j.leary/py_envs/scLANE_env2
 version:        3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
 numpy:          /home/j.leary/.local/lib/python3.10/site-packages/numpy
 numpy_version:  1.24.2
 
 NOTE: Python version was forced by use_python() function

──────────────────────────────────────────────────────────────────────────────